clear all;
clc;

xe = 0;
ye = 3;
ze = 0;

cd('patient1_8_12_13');
xyz_dc = load('p1_dc_xyzloc.txt');
xyz_dr = load('p1_dr_xyzloc.txt');


[noaxon,junk] = size(xyz_dc);

r_dc = sqrt((xyz_dc(:,1)-xe).^2+(xyz_dc(:,2)-ye).^2+(xyz_dc(:,3)-ze).^2);
r_dr = sqrt((xyz_dr(:,1)-xe).^2+(xyz_dr(:,2)-ye).^2+(xyz_dr(:,3)-ze).^2);

[rr_dc,si_dc] = sort(r_dc);
[rr_dr,si_dr] = sort(r_dr);

% debugging 1
% figure;
% plot(xyz_dc(:,1),xyz_dc(:,2),'k.');
% hold on;
% for ii = 1:noaxon
%     plot(xyz_dc(si_dc(ii),1),xyz_dc(si_dc(ii),2),'ro');
%     gtmp = input('');
% end
% hold off;

% ------------------------- Importing Data --------------------------------

cd('transtripole_1ies');
tt1phi_dc = load('sub_TT1_dc_phi_p1_d0p2_t0.txt');
tt1phi_dc = tt1phi_dc(:,1:8:end);
tt1phi_dc = tt1phi_dc(si_dc,:);
tt1sd_dc = diff(tt1phi_dc,2,2);
tt1msd_dc = max(tt1sd_dc,[],2);

tt1phi_dr = load('sub_TT1_dr_phi_p1_d0p2_t0.txt');
tt1phi_dr = tt1phi_dr(:,1:8:end);
tt1phi_dr = tt1phi_dr(si_dr,:);
tt1sd_dr = diff(tt1phi_dr,2,2);
tt1msd_dr = max(tt1sd_dr,[],2);
cd ..;

cd('transtripole_3ies');
tt3phi_dc = load('sub_TT3_dc_phi_p1_d0p2_t0.txt');
tt3phi_dc = tt3phi_dc(:,1:8:end);
tt3phi_dc = tt3phi_dc(si_dc,:);
tt3sd_dc = diff(tt3phi_dc,2,2);
tt3msd_dc = max(tt3sd_dc,[],2);

tt3phi_dr = load('sub_TT3_dr_phi_p1_d0p2_t0.txt');
tt3phi_dr = tt3phi_dr(:,1:8:end);
tt3phi_dr = tt3phi_dr(si_dr,:);
tt3sd_dr = diff(tt3phi_dr,2,2);
tt3msd_dr = max(tt3sd_dr,[],2);
cd ..;

cd('longtripole_1p5ies');
lt1p5phi_dc = load('sub_LT1p5_dc_phi_p1_d0p2_t0.txt');
lt1p5phi_dc = lt1p5phi_dc(:,1:8:end);
lt1p5phi_dc = lt1p5phi_dc(si_dc,:);
lt1p5sd_dc = diff(lt1p5phi_dc,2,2);
lt1p5msd_dc = max(lt1p5sd_dc,[],2);

lt1p5phi_dr = load('sub_LT1p5_dr_phi_p1_d0p2_t0.txt');
lt1p5phi_dr = lt1p5phi_dr(:,1:8:end);
lt1p5phi_dr = lt1p5phi_dr(si_dc,:);
lt1p5sd_dr = diff(lt1p5phi_dr,2,2);
lt1p5msd_dr = max(lt1p5sd_dr,[],2);
cd ..;

cd('longtripole_6ies');
lt6phi_dc = load('sub_LT6_dc_phi_p1_d0p2_t0.txt');
lt6phi_dc = lt6phi_dc(:,1:8:end);
lt6phi_dc = lt6phi_dc(si_dc,:);
lt6sd_dc = diff(lt6phi_dc,2,2);
lt6msd_dc = max(lt6sd_dc,[],2);

lt6phi_dr = load('sub_LT6_dr_phi_p1_d0p2_t0.txt');
lt6phi_dr = lt6phi_dr(:,1:8:end);
lt6phi_dr = lt6phi_dr(si_dc,:);
lt6sd_dr = diff(lt6phi_dr,2,2);
lt6msd_dr = max(lt6sd_dr,[],2);
cd ..;

cd('other_geom');
atpphi_dc = load('sub_atp_dc_phi_p1_d0p2_t0.txt');
atpphi_dc = atpphi_dc(:,1:8:end);
atpphi_dc = atpphi_dc(si_dc,:);
atpsd_dc = diff(atpphi_dc,2,2);
atpmsd_dc = max(atpsd_dc,[],2);

atpphi_dr = load('sub_atp_dr_phi_p1_d0p2_t0.txt');
atpphi_dr = atpphi_dr(:,1:8:end);
atpphi_dr = atpphi_dr(si_dc,:);
atpsd_dr = diff(atpphi_dr,2,2);
atpmsd_dr = max(atpsd_dr,[],2);
cd ..;

cd ..;

% ------------------------- DC/DR Rangles ---------------------------------

lt1p5_dc_range=[min(lt1p5msd_dc),max(lt1p5msd_dc)];
lt1p5_dr_range=[min(lt1p5msd_dr),max(lt1p5msd_dr)];

lt6_dc_range=[min(lt6msd_dc),max(lt6msd_dc)];
lt6_dr_range=[min(lt6msd_dr),max(lt6msd_dr)];

tt1_dc_range=[min(tt1msd_dc),max(tt1msd_dc)];
tt1_dr_range=[min(tt1msd_dr),max(tt1msd_dr)];

at_dc_range=[min(atpmsd_dc),max(atpmsd_dc)];
at_dr_range=[min(atpmsd_dr),max(atpmsd_dr)];

figure;
hold on;
plot([0.95,0.95],lt6_dc_range,'k-','LineWidth',4);
plot([1.05,1.05],lt6_dr_range,'k--','LineWidth',4);
plot([1.95,1.95],tt1_dc_range,'k-','LineWidth',4);
plot([2.05,2.05],tt1_dr_range,'k--','LineWidth',4);
plot([2.95,2.95],at_dc_range,'k-','LineWidth',4);
plot([3.05,3.05],at_dr_range,'k--','LineWidth',4);
hold off;
ylabel('Maximum \delta^2\Phi (V)','FontSize',30);
set(gca,'FontSize',30,'XTick',1:1:3);
ylim([0,0.75]);
xlim([0.75,3.25]);

% ------------------------- Additional Processing -------------------------

allsd_dr = cat(2,tt1msd_dr,lt1p5msd_dr,atpmsd_dr);

i50 = noaxon/2;
maxsd50_dc = [tt1msd_dc(i50),lt1p5msd_dc(i50),atpmsd_dc(i50)];

% ------------------------- Decay of SD -----------------------------------

figure;

ss=1:200;%[1,10:10:200];

% Maximum SD across DC
yo = 0.001;
yf = 0.7;
subplot(2,2,[1,3]);
hold on;
plot(100*ss/max(ss),-sort(-tt1msd_dc),'r-','LineWidth',4);
plot(100*ss/max(ss),-sort(-lt1p5msd_dc)...
    ,'b-','LineWidth',3,'MarkerSize',30);
hold off;
ylim([yo,yf]);
xlabel('% DC Fibers','FontSize',30);
ylabel('Maximum \delta^2\Phi (V)','FontSize',30);
set(gca,'FontSize',30,'YScale','Linear','YTick',0:0.1:yf,'XTick',[10,25,50,75,90,100]);
legend('TT, IES = 1 mm','LT, IES = 1.5 mm','E');

Q = tt1msd_dr;
R = lt1p5msd_dr;

subplot(2,2,2);
hold on;
[a,b]=hist(Q,40);
bar(b,100*a/sum(a),'BarWidth',1);
ten_1=sort(tt1msd_dc);
ten_1=ten_1(20);
plot([ten_1,ten_1],[0,55],'k--','LineWidth',4);
hold off;
ylabel('% DR Fibers','FontSize',30);
h = findobj(gca,'Type','patch');
set(h,'FaceColor','r','EdgeColor','w','facealpha',0.75);
set(gca,'FontSize',30,'YTick',[10,25,50]);
xlim([0,0.17]);
ylim([0,55]);

subplot(2,2,4);
hold on;
aa=hist(R,b);
bar(b,100*aa/sum(aa),'BarWidth',1);
ten_2=sort(lt1p5msd_dr);
ten_2=ten_2(20);
plot([ten_2,ten_2],[0,15],'k--','LineWidth',4);
hold off;
xlabel('Maximum \delta^2\Phi (V)','FontSize',30);
ylabel('% DR Fibers','FontSize',30);
h1 = findobj(gca,'Type','patch');
set(h1,'FaceColor','b','EdgeColor','w','facealpha',0.75);
set(gca,'FontSize',30,'YTick',[10,15]);
xlim([0,0.17]);
ylim([0,15]);

% ------------------------- Plotting Results ------------------------------



return

% Maximum SD across DR
gm1='TT, IES = 1 mm';
gm2='LT, IES = 1.5 mm';
gm3='ATP';
geomlab = {gm1,gm2,gm3};

figure;
hold on;
h2 = boxplot(allsd_dr,'labels',geomlab,'Color','k');
plot(1:3,maxsd50_dc,'kd','MarkerSize',15,'MarkerFaceColor','k');
title('Dorsal Root (DR) Fibers','FontSize',30);
ylabel('Maximum \Delta^2\Phi (V)','FontSize',30);
legend('DC_{50}','Location','NE');
set(gca,'FontSize',24);
set(gca,'xtick',1:size(allsd_dr,2),'xticklabel',geomlab,'FontSize',26);
set(h2,'LineWidth',3);


